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Abstract 

In a plethora of applications dealing with inverse problems, e.g. in image processing, social net¬ 
works, compressive sensing, biological data processing etc., the signal of interest is known to be 
structured in several ways at the same time. This premise has recently guided the research to the 
innovative and meaningful idea of imposing multiple constraints on the unknown parameters involved 
in the problem under study. For instance, when dealing with problems whose unknown parameters form 
sparse and low-rank matrices, the adoption of suitably combined constraints imposing sparsity and low¬ 
rankness, is expected to yield substantially enhanced estimation results. In this paper, we address the 
spectral unmixing problem in hyperspectral images. Specifically, two novel unmixing algorithms are 
introduced, in an attempt to exploit both spatial correlation and sparse representation of pixels lying in 
homogeneous regions of hyperspectral images. To this end, a novel mixed penalty term is first defined 
consisting of the sum of the weighted and the weighted nuclear norm of the abundance matrix 
corresponding to a small area of the image determined by a sliding square window. This penalty term 
is then used to regularize a conventional quadratic cost function and impose simultaneously sparsity 
and row-rankness on the abundance matrix. The resulting regularized cost function is minimized by a) 
an incremental proximal sparse and low-rank unmixing algorithm and b) an algorithm based on the 
alternating minimization method of multipliers (ADMM). The effectiveness of the proposed algorithms 
is illustrated in experiments conducted both on simulated and real data. 

P. V. Giampouras, K. E. Themelis, A. A. Rontogiannis and K. D. Koutroumbas are with the Institute for Astronomy, 
Astrophysics, Space Applications and Remote Sensing (lAASARS), National Observatory of Athens, I. Metaxa & Vas. Pavlou 
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Index Terms 


Semi-supervised spectral unmixing, hyperspectral images, simultaneously sparse and low-rank ma¬ 
trices, proximal methods, alternating direction method of multipliers (ADMM), abundance estimation 


1. Introduction 

Spectral unmixing (SU) of hyperspectral images (HSIs) has attracted considerable attention 
in recent years both in research and applications. SU can be considered as the process of a) 
identifying the spectral signatures of the materials (endmembers) whose mixing generates the 
(so called) mixed pixels of an HSI and b) deriving their corresponding fractions {abundances) in 
the formation of each HSI pixel, Q. The latter constitute the so called abundance vector of the 
pixel. This two step procedure has given rise to a plethora of methods tackling either one or both 
these two tasks. Diverse statistical and geometrical approaches have been lately put forward in 
literature addressing the first step, commonly known as endmembers’ extraction (e.g. 0 , 0 ). 
On the other hand, there have been many research works that assume that the spectral signatures 
of the endmembers are available and focus on the abundance estimation task. Algorithms that fall 
into this class, need to make a fundamental assumption concerning the inherent mixing process 
that generates the spectral signatures of the HSI pixels. 

In view of the latter, the linear mixing model (LMM) holds a dominant position being widely 
adopted in numerous state-of-the-art unmixing algorithms (see e.g. O and the references therein). 
More specifically, these algorithms are based on the premise that the pixels’ spectral signatures 
are generated by a linear combination of endmembers’ spectra contained in a predefined set, 
usually termed as endmembers’ dictionary. Abundance estimation is henceforth treated as a linear 
regression problem. LMM has prevailed over other models, due to its simplicity and mathemat¬ 
ical tractability. Physical considerations that naturally arise, impose various constraints on the 
unmixing problem. In this context, the so-called abundance nonnegativity and the abundance 
sum-to-one constraints are usually adopted. That said, unmixing can be viewed as a constrained 
linear regression problem. 

In an attempt to achieve better abundance estimation results, recent novel ideas promote the 
incorporation of further prior knowledge in the unmixing problem. In light of this, several 
methods bring into play the sparsity assumption, 0 - 0 . Its adoption is justified by the fact that 
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only a few of the available endmembers participate in the formation of a given pixel, especially 
in the case of large size endmembers’ dictionaries. Put it in other terms, it is envisaged that 
pixels’ spectral signatures accept sparse representations with respect to a given endmembers’ 
dictionary. Furthermore, one could also say that the abundance vectors corresponding to the 
pixels of HSIs are deemed having only a few non-zero values. Practically speaking, sparsity 
is imposed on abundances by means of norm regularization, 0-0 when a deterministic 
approach is followed. On the other hand, in Bayesian schemes appropriate sparsity inducing 
priors are adopted for the abundance vectors, Q, 0. Spatial correlation is another constraint 
that has recently been incorporated in the unmixing process, offering stimulating results, Q-[ 111. 
In that vein, the additional information that exists in homogeneous regions of HSIs is subject to 
exploitation. Actually, in such regions, there is a high degree of correlation among the spectral 
signatures of neighboring pixels. It is hence anticipated that there should also be correlation 
among the abundance vectors corresponding to these pixels. This has led to the development of 
novel unmixing schemes, whereby the information provided by the neighboring pixels is taken 
into account in the abundance estimation of each single pixel. 

In this spirit, a collaborative deterministic scheme, termed CLSUnSAL, was recently proposed 
in GD. which uses a wealth of information stemming from all the pixels of the examined 
HSI. CLSUnSAL adopts dictionaries consisting of a large amount of endmembers. Then it 
assumes that spatial correlation translates into abundance vectors sharing the same support set i.e., 
presenting a similar sparsity pattern. Thus, the matrix whose columns are the abundance vectors 
of all HSI pixels (called abundance matrix) should meaningfully be of a joint-sparse structure 
To impose joint-sparsity, CLSUnSAL applies a £ 2,1 norm on the sought abundance matrix, which 
is then used to penalize a suitably defined quadratic cost function. Minimization of the resulting 
regularized cost function is performed by an alternating direction method of multipliers (ADMM), 
0. A similar perspective is followed in p0|, however in a “localized” fashion. Specifically, 


1101 proposes the use of a 3 x 3 square window that slides all over the image. The abundance 
vector of the central pixel is then inferred, by taking into account the spectral signatures of the 
adjacent pixels contained in the window. Based on this idea, two algorithms are derived. First the 
MMV-ADMM, which in a similar to CLSUnSAL fashion, seeks joint-sparse abundance matrices 


'a joint-sparse Bayesian unmixing scheme has been also presented in jizj. 
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Utilizing the £ 2,1 norm, and second the LRR algorithm that promotes a low-rank structure on 
the abundance matrix. Actually, the LRR algorithm presents an alternative way of modelling 
the spatial correlation among neighboring pixels. That is, it assumes that the correlation among 
pixels’ spectral signatures is reflected as linear dependence among their corresponding abundance 
vectors. Apparently, the matrix formed by these abundance vectors should be of low rank. That 
said, a nuclear norm is imposed on the abundance matrix, and a properly adapted augmented 
Lagrangian cost function is minimized in an alternating minimization fashion. 

In this paper, we introduce a novel idea for performing abundance estimation in HSIs under 
the LMM, that simultaneously takes spatial correlation and sparsity into consideration. Similarly 


to [101, we utili z e a. k x k square sliding window with k odd, and we consider the spectral 
signatures of adjacent pixels lying in it. Departing from the usual paradigm, we propose to 
seek for K^-column abundance matrices that are simultaneously sparse and low-rank. SU is thus 


formulated as a sparse reduced-rank regression problem, [14|. As stated earlier, low-rankness 
arises naturally in abundance matrices corresponding to relatively homogeneous regions, due 
to the linear dependence of the respective abundance vectors. At the same time, sparsity is a 
reasonable hypothesis that still holds independently, as explained above, within each individual 
abundance vector. Broadly speaking, imposing multiple structures on the same mathematical ob¬ 
ject when dealing with inverse problems is a strategy still in its very infancy in signal processing 


and machine learning literature, [ 15|-[ 18|. The aforementioned sparsity and low-rank constraints 
give rise to a mixed penalty term that regularizes a least squares fitting function through the 
weighted £1 norm and the weighted trace norm of the abundance matrices, respectively. In order to 
minimize the cost function, two novel iterative algorithms are proposed, namely the incremental 


proximal sparse low-rank unmixing algorithm (IPSpLRU), inspired by [ [T9| , and the alternating 
direction sparse and low-rank unmixing algorithm (ADSpLRU). As implied by their names, 
IPSpLRU is based on proximal operators of the individual terms that compose the cost function. 
On the other hand, ADSpLRU is an ADMM based approach properly adapted to our problem 
formulation. The proposed algorithms are compared with state-of-the-art unmixing techniques 
and their effectiveness is demonstrated via extensive simulated and real-data experiments. 

Notation: Matrices are represented as boldface uppercase letters, e.g., X, and, column vectors 
as boldface lowercase letters, e.g., x, while the i-th component of vector x is denoted by Xi 
and the ij-th element of matrix X by Xij. Moreover, ^ denotes transposition, I at is the N x N 
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PIXELS SIGNATURES Y ENDMEMBERS DICTIONARY (p 





Fig. 1: Graphical illustration of the sliding window approach of our unmixing algorithms. 
Abundance matrix is considered sparse and low-rank (in this example rank = 2). Blue cells 
in matrix W represent zero values. 


identity matrix and 0 is a zero matrix with respective dimensions, 1 denotes the all ones vector, 
rank(X) is the rank of X, tr[X] denotes the trace of matrix X, diag(x) is a diagonal matrix 
with the elements of vector x on its diagonal, aj(X) is the i\h largest singular value of X, 11-112 
is the standard £2 (Euclidean) vector norm, ||X||* = Tr(VX'^X) = denotes 

the nuclear norm (or trace norm), ||X||i = *^he sum of the absolute values of 

all entries of X (called the £i norm), ||X||ir = stands for the Frobenius norm. 

AA(-) denotes the Gaussian distribution. Also, stands for the fc-dimensional Euclidean space 
and denotes the fc-dimensional non-negative orthant. The matrix inequality X > Y declares 
element-wise operation and 0 stands for component-wise multiplication between matrices of the 
same size. 


IL Problem Formulation 

We consider an L-spectral band hyperspectral image, with each of its pixels being composed 
of N endmembers. Eet $ = <^ 2 , ■ ■ ■, <Pn] stand for the Lx N endmembers’ dictionary, 

where cp^ G M+, i = 1,2,..., N, is the spectral signature of the Ah endmember. Consider also a 
small sliding square window that contains K adjacent pixels (K = kx n), with the measurement 
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spectra y^, k = 1,2,... ,K, that are assumed to share the same endmember matrix as shown 
graphically in Fig. for TF = 9. In matrix notation, let Y = [yi,y 2 , • • • ,yir] be the L x iF 
matrix containing the spectra of the K pixels in the window as its columns. Utilizing the linear 
mixing model (LMM), the mixing process can be described by the equation 


Y = $W + E, 


( 1 ) 


where W e is the abundance matrix whose columns are the iV-dimensional abundance 

vectors of the corresponding K pixels, and E G is an i.i.d., zero-mean Gaussian noise 

matrix. Due to physical considerations, the abundance coefficients in W should satisfy two 
constraints, namely, the abundance nonnegativity and the abundance sum-to-one constraints. 


|20|, i.e.. 


W > 0, and l^W = 1^. 


( 2 ) 


Nevertheless, in the following we relax the sum-to-one constraint based on the reasoning pre¬ 
sented in Q. That said, the general problem considered in this paper is the following: “given 
the spectral measurements Y and the endmember matrix $, estimate the abundance matrix 
W subject to the nonnegativity constraint”. This is a typical inverse problem, which has been 
addressed via many methods in the signal processing literature. However, the efficacy of the 
proposed approach lies on the exploitation of the intrinsic structural characteristics of W, i.e., 
sparsity and low-rankness. To this end, we impose concurrently two naturally justified structural 
constraints on the abundance matrix W, that promote low-rankness and sparsity. 

Low-rankness property: A logical consideration is that all pixels belonging to the same window 
are correlated, i.e., they are composed of the same materials, although maybe in different propor¬ 
tions. This property suggests that the abundance matrix W to be estimated has linearly dependent 
columns and thus is either low-rank, or it can be well-approximated by a low-rank matrix. 
In the bibliography, low-rank matrix estimation techniques have recently emerged as powerful 
estimation tools, e.g., p0|, pT|-[|23|. These estimators are mainly based on regularization by 


the nuclear norm of W. A similar regularization is also adopted in this paper in order to impose 
the low-rank constraint. 

Sparsity property: Another typical assumption is that only a small portion of the N endmem- 
bers will be present in the spatial area marked by the n x n shifting window. In other words. 
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it is safe to assume that the abundanee matrix W has a sparse representation in terms of the 
endmember matrix dietionary This motivates the use of a sparsity-eognizant estimator for the 
abundanee matrix W, whieh is envisaged to produee more robust unmixing results. It should be 
noted that sparsity has already been sueeessfully exploited in many speetral unmixing algorithms, 

e.g., 0-®, [|TT|, 

It is worth mentioning that the sparsity of W does by no means invalidate its low-rankness. 
On the contrary, both structural hypotheses on W are assumed to hold simultaneously, although 
low-rankness implicitly imposes some kind of structure on sparsity. So far, reports in the spectral 
unmixing literature explore either the sparsity, e.g. Q. d). or the low-rankness property of W, 


e.g. [ 10|. This is the first time, to the best of our knowledge, that spectral unmixing is formulated 
as a simultaneously sparse and low-rank matrix estimation problem. That is, we seek a matrix 
W > 0 that, apart from fitting the data well in the least squares sense, it has minimum rank and 
only a few positive elements. To achieve this, we define the following optimization problem. 


( 3 ) 


(PI): W=argmin -||Y - $W||^ + 7 ||W||i + r||W| 


where 7 , r > 0 are parameters that control the trade-off between the sparsity and rank regular¬ 
ization terms and the data fidelity term. Being parametrized, (PI) becomes flexible enough to 
impose either one of the two structures on W. For example, by setting 7 = 0, (PI) results in 
searching for a matrix that is of low-rank structure. Accordingly, setting r = 0 is tantamount 
to searching for a sparse matrix. The flexibility of the proposed model provides certainly an 
advantage over either low-rank or sparse estimation methods, as it will also be demonstrated 
later, in the experimental results section. 

It is also worth pointing out that (PI) involves the convex surrogates of the zero norm ||W||o 
and rank(W), i.e., the £1 and the nuclear norm, respectively. In an attempt to promote further 
the robustness and consistency of the proposed estimator, we propose to use weighted l\ and 
nuclear norms in (PI). Such an approach is expected to enhance the sparsity on the individual 


elements and the singular values (Ti(W), e.g. [25 |-p^- These weighted norms are defined 
as 


N K 


|A 0 Will = y^^y^^aij\wij\, 

i=l j=l 


( 4 ) 
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where and bi are nonnegative weighting eoeffieients (i.e. > 0 and > 0). Utilizing Q 

and Q, the proposed optimization problem is rewritten as 

(P2) ; W=argmin (i|| Y - $W||| + 7|| A © W||i + r|| W||b,4 • (6) 

To the best of our knowledge, sueh a formulation, has not been used before as a regularizer 
for promoting simultaneously sparsity and low-rankness. In the following, problem (P2) will be 
studied under the following assumption. 

Assumption 1: For the weighting eoeffieients bi of the nuelear norm it holds that bi = b. 




Under Assumption 1, the nuelear norm is eonvex [271, [29 j, while the weighted ii norm is always 
eonvex for nonnegative A. Thus, the overall cost function of (P2) is convex. Alternative options 
for the selection of parameters A and b are discussed in section III-C. Although convex, (P2) 
under Assumption 1 is a nontrivial problem to solve, due to the non-differentiable form of the ii 
and nuclear norm regularizers, [ [30| . In the following, we suitably explore two standard convex 
optimization tools to tackle this problem; an incremental proximal method and an alternating 
direction method of multipliers (ADMM) based technique. 


III. The Proposed algorithms 

In this section we present two algorithms to address the non-smooth, constrained, convex 
optimization problem in (P2). The first one is an incremental proximal algorithm, recently 


presented and analyzed in [19[, which makes use of the proximal operators of all the terms 
appearing in (P2), while the second exploits the splitting strategy of the ADMM philosophy. 


[13|. 


A. Incremental proximal descent sparse and low-rank unmixing algorithm 

Let us first recall that the proximal operator of a function /(■) is defined as [31 [, [32[, 

P’^o^A/{.)(^) = argmm |^/(x) + - d||2^ , ( 7 ) 
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where v and x G dom(/), the domain of /. In |19| the following minimization problem 
is eonsidered 


mm 

xeA’ 




( 8 ) 


i=l 


where /i(x),i = are eonvex funetions and A” C is a elosed eonvex set. One 


version of the algorithm proposed in [19| to solve this problem is the following. The proximal 
operators of all //s are first derived and then a sequential seheme is defined, in whieh the 
proximal operator of /i(x) is evaluated at the point provided by its predeeessor (the proximal 
operator of /j_i(x)), for i = 2,3,..., m. This proeedure is repeated in a eyelie manner at eaeh 
iteration of the algorithnj^ A eonvergenee and rate of eonvergenee analysis of this ineremental 


proximal seheme is also given in [ 19|. After this short introduetion, we may observe that (P2) in 
(|^ with b = 61 has exaetly the same form with the minimization problem in Q, with respeet to 
W. Embedding the nonnegativity to the eost funetion in (|^ we obtain the following regularized 
quadratie loss funetion, 


A(W) = -||Y 


$W||| 


|A©W||i + r||W||b.*+XK4W), 


(9) 


where the nonnegativity eonstraint is now replaeed by the (eonvex) indieator funetion X]r_^(W), 
whieh is zero when all > 0, i = 1,2,..., Y, j = 1, 2,..., K, and +oo if at least one Wij is 
negative. Typieally, we wish to minimize Xi(W) with respeet to W. Notiee that Xi(W) is the 


sum of four eonvex funetions and the ineremental proximal algorithm of [19| ean be applied 
direetly in our problem. Next, the proximal operators of all four eonvex funetions are obtained. 
Starting with the least squares fitting term, we readily get 


prox;^i|IY„$.||2 (W) = ($^ $ + Y + A-^W) 


( 10 ) 


Before we give the proximal operators for the next three terms, some neeessary definitions are 
in order. First, we define the soft-thresholding operator on matrix W = [w]ij as 

SHRa(W) = sign(W) max(0, |W| - A), (11) 

where A = [5] is the matrix that eontains thresholding parameters. Note that the soft-thresholding 


in (111 is performed in an element-wise manner, i.e., SHR 5 .^,(t(;jj) = sign{wij) max(0, \xij\—6ij). 

^Instead of sequential, a randomized evaluation of the proximals of fi’s could be also employed, jl9|. 
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Notably, when we apply the soft-thresholding operator on a diagonal matrix, we shrink only the 
elements belonging to its diagonal. These elements are assumed to be shrinked by thresholding 
parameters eontained in a veetor. With this in mind, we define the singular value thresholding 
operation by 

SVTs{W) = USHR5(S) 

where W = USV^ is the singular value deeomposition (SVD) of W, and S is the veetor 
whose entries are the thresholding parameters that reduee the eorresponding diagonal elements 
of matrix S. Finally, we define the projeetion operator on the set of nonnegative real numbers, 

f 0, V < 0 

nR, (n) = argminlx — = < , (12) 

xeR+ ^ > 0 

whieh ean also be applied to matriees in an element-wise manner. 

Utilizing the above definitions, we ean eompute the proximal operators for all regularizing 
eonvex funetions in (|^. Speeifieally, prox^|j^Q.||^(W) is eomputed by soft-thresholding matrix 
W with yA as follows. 


prox.^ljA0.|l^(W) = SHR.,a(W). (13) 

Similarly, the proximal operator of the nuelear norm ean be expressed via a soft thresholding 
operation on the singular values of W, i.e.. 


prox,||.|l^^(W) = SVT,b(W). (14) 

Moreover, the eomputation of proxjj|j^(.)(W) reduees to a projeetion operation, i.e., 

proxij^^(.)(W) = nR^(W). (15) 

The proposed ineremental proximal sparse and low-rank unmixing algorithm (IPSpLRU) iterates 
among the proximal operators ( [T^ , ( [13] ), ( [T4| ) and ( [T5] ) in a eyelie order until eonvergenee, p9|. 
IPSpLRU is summarized in Algorithm Note that in order to retain the eonvexity of the 
eomposite funetions, the weighting parameters A and b are initialized and kept fixed during the 
exeeution of the algorithm. The issue of dynamie seleetion of these parameters is diseussed in 


Seetion |III-C[ Coneeming the eomputational eomplexity of IPSpLRU, the most eomplex step 
is the SVD of the abundanee matrix W*, whieh takes plaee at eaeh iteration and is of the order 
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Algorithm 1 The proposed IPSpLRU algorithm 
Inputs Y, $ 

Seleet parameters A, b, A, r, 7 

Set t = 0, R = P = ^^Y, Q = RP 

Initialize W° 

repeat 

W* = Q + A-^RW* 

W* = prox,||AQ.||^(W‘) 

W* = prox^||.||^^(W*) 

= proxj^^(.)(W*) 

until eonvergenee 

Output : Abundanee matrix W = 


of 0{KN^ + K^), [331. Note that matriees R = + X-^1n)~\ P = $^Y and Q = RP 


are eomputed only onee at the initialization stage and thus the first step of the algorithm just 
requires a fast matrix-by-matrix multiplieation. The algorithm eonverges rapidly and terminates 
when either the following stopping eriterion is satisfied, 


ll^m _ w*||2, 
l|W*||^ 


< 6 


(16) 


where 5 is a predefined threshold value, or a preset maximum number of iterations is reaehed. In 
the following seetion we present an alternative approaeh to solve the same problem by employing 
a primal-dual ADMM type teehnique. 


B. Alternating direction method of multipliers for sparse and low rank unmixing 

In this seetion, we develop an instanee of the alternating direetion method of multipliers that 
solves the abundanee matrix estimation problem (P2). To proeeed, we utilize the auxiliary matrix 
variables 171 , 172,^23 and 174 of proper dimensions (similar to [[TT|, p^), and reformulate the 


original problem (P2) into its equivalent ADMM form, [13|, 
(P3): 


mm < -||17i 
,r22,r23,r24 [ 2 


YI |2 
1 IIf 


IIA © 172111 + 'r||173||b,* + 7*^( 174 ) ^ (17) 

s.t. 17i — $W = 0, 172 — W = 0, 173 — W = 0, 174 — W = 0 
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Based on (P3), the following augmented Lagrangian funetion is defined, 


-^2(W, fii, f 22 , fis; ^^4) — 2 ~ + 7II A 0 f22|| 1 + 'rllfisllb,* + 

+ tr [Df {Cli - $W)] + tr [D^ (f22 - W)] + tr [D[ (fig - W)] + tr [D[ (fl4 - W 


^ (||$W - flilll + ||W - fl 2 ||?. + ||W - flgll^ + ||W - fl4||?.) 


(18) 


where the L x K matrix Di, and the N x K matriees D 2 , Dg, D 4 are the Lagrange multipliers 
and /i > 0 is a positive penalty parameter. Note that again the nonnegative weights A and b are 
eonsidered to be eonstant and assumption 1 also holds here. Let 


n = 


’ f7i ' 




’ -II 

0 

0 

0 

172 

, G = 

liv 

, and B = 

0 

— '^N 

0 

0 

17g 


Iat 


0 

0 

—In 

0 

174 


liv 


0 

0 

0 

—Ixf 


Then (18) ean be written in an equivalent form as 


(19) 


1^3(W, fl, A) — 2 11^1 ~ + 7 IIA O r22||i + I'llflslib,♦ 

+ XM^(n4) + |||GW + Bfl - A|||, 


( 20 ) 


where A = [Af A^ Ag Aj]^, A* = (l//i)Di,z = 1,... ,4, eontains the sealed Lagrange mul¬ 


tipliers. Having expressed the augmented Lagrangian funetion as in (20), the ADMM proeeeds 
by minimizing £g(W, r 2 ,A) sequentially, eaeh time with respeet to a single matrix variable, 
keeping the remaining variables at their latest values. The dual variables (Lagrange multipliers) 
are also updated via a gradient aseend step at the end of eaeh alternating minimization eyele. 

To elaborate further on the steps of the ADMM, the optimization with respeet to W gives 


= argrmn £g(W, f 2 *, A*) 

= + 3Ijv)”' + \\) + nl + K\ + il\ + K\ + n\ + A*] 


( 21 ) 


Next, the optimization with respeet to fii is performed as 

- o T.gv.n 1 v-i ^ 

1 


= argmin £g(W*^^, 17, A* 


1 + /i 


(Y + /i($W*+i - A*)) . 


( 22 ) 
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The remaining auxiliary variables f22, ^ 3 , and are involved in non-differentiable norms, 
namely, the weighted ii norm, the weighted nuelear norm, and the indicator function, respec¬ 
tively. In this regard, the minimization task with respect to these variables resolves to computing 
some of the proximity operators that we introduced in the previous section. Minimizing ( |20| ) 
with respect to 0,2 yields 

= argmin A*) 

= SHR.^a(W*+^ - A*). (23) 

In the same vein, 0.3 is computed by a shrinkage operation, 

= argmin £ 3 (W*’''^, f2, A*) 

fiS 

= SVT,b(W^+^ - A*). (24) 


Next, for the auxiliary variable f 24 , a projection onto the nonnegative orthant is required, 

= argmin £ 3 (W‘+\ 12, A*) 




= nK^(W*+'-A‘). (25) 

At the final step of the proposed method, the scaled Lagrange multipliers in A are sequentially 


updated by performing gradient ascent on the dual problem pA| |, as follows, 

A*+^ = A* - 

A*+^ = A* - + Ol+\ i = 2,3,4 


(26) 


The proposed algorithm, termed the alternating direction sparse and low-rank unmixing algorithm 
(ADSpLRU), is summarized in Algorithm An iteration of ADSpLRU consists of the update 
steps given in ( [^ , ( |^ , ( [^ , ( |^ , and ( [^ . Its computational complexity is 0{KLN + 


KN"^) per iteration, slightly higher than that of IPSpLRU, since it usually holds L > N. However, 
as verified by the simulations of the next section, ADSpLRU converges a little faster than 
IPSpLRU to a slightly lower steady-state erroij^ while its convergence is also guaranteed as 
explained in [34|. Note that all functions that form the objective function £i(W) in Q, are 


^The reason for this may be that ADSpLRU manipulates the whole cost function at each step, while IPSpLRU splits the cost 
function in a number of convex terms and treats each term individually at every step of the algorithm. 
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Algorithm 2 The proposed ADSpLRU algorithm 
Inputs Y, $ 

Seleet parameters A, b, /x, r, 7 
Set t = 0, R = + 3Ijv)“^ 

Initialize W°, 12°, A° 

repeat 

W*+i = + A*) + 12^ + K\ + il\ + K\ + n\ + A*] 

12*+^ = 1/(1 + /i) (Y + /i - AO) 

12*+^ = SHR^a(W*+i - A*) 

120^ = SVT,b(W*+i - A*) 

120 '= Rk^W^+i - AO 
AO' = A* - $W*+' + I2O' 

AO' = A* - + I 2 O', i = 2,3,4 

until eonvergenee 
Output : Abundanee matrix W = 


elosed, proper and eonvex. Sinee matrix G has full eolumn rank, the eonvergenee eonditions 


defined in [34| are met and if an optimal solution exists, ADSpLRU eonverges, for any /i > 0. 
This in turn implies that for the primal and dual residuals r*, d* given by 


r* = GW' + B12', 

d' = fiG^B ( 12 ' - 12 '-') 

it holds that, r' —)■ 0 and d' —0, respeetively, as t ^ 00 . In this work, ADSpLRU stops when 
either the following termination eriterion 


lir'ils < C and ||d '||2 < C (27) 

holds for the primal and dual residuals, where ( = 7 / {3N + |j^, (the relative toleranee 

Qrei ^ g value depending on the applieation, and in our experimental study has been 

empirieally determined to 10 “^), or the maximum number of iterations is reaehed. 
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C. Selection of weighting coefficients and regularization parameters 

As mentioned previously, in both IPSpLRU and ADSpLRU the weighting coefficients A and 
b are predetermined, satisfy certain constraints and remain constant during the execution of the 


algorithms. As is widely known, p5| , p 6 | , p9| |, a proper selection of these parameters is quite 
crucial as for the accuracy of the estimations. In view of this, two potential choices are a) to 
select the weighting coefficients based on the least squares estimate of W i.e., 

1 \ / 1 


Oij — 


,LS 


w-f + e 


and bj = 


(28) 


a,(W^^) + e, 

where e = 10 “^® is a small constant added to avoid singularities and b) to update them at each 
iteration t of the algorithms based on the current estimate W* of W i.e., 

1 \ / 1 




wjj + e 


and 6 ; = 


(29) 


It should be noted, that both these two options render the minimization problem (P2) non- 
convex, since the weighted nuclear norm is known to be convex only if the weights bi, i = 


1,2,... ,K are nonnegative and non-ascending, [29|, [351. Additionally, the reweighting norms 


minimization problem is known to be inherently nonconvex, [26|, while its theoretical conver¬ 
gence analysis for these cases is difficult to be established. Nevertheless, numerous research 
works advocate the positive impact of these nonconvex weighted norms on the performance of 


constrained estimation tasks [26|, [27|, [29|, [351. Along this line of thought, the algorithms 
presented in the previous section are modified by adopting the reweighting scheme given by 
( [29| ). As verified in our empirical study presented in the next section, such an option enhances 
to a large degree the effectiveness of the proposed algorithms, while no numerical issues have 
been encountered in our experiments. 

As far as the remaining parameters is concerned, A and /i, which control the convergence 
behavior of IPSpLRU and ADSpLRU, respectively, take positive values, with p close to zero 
and A on the order of 1. In all our experiments we fixed p = 0.01 and A = 0.5. On the other 
hand, the low-rank and sparsity promoting parameters r and 7 are chosen via fine-tuning, as is 
commonly done in relevant deterministic schemes. This is so because the optimal set of these 
parameters depends on the unknown in advance particular structure of the sought abundance 
matrix, an issue which is further explained in the next section. 
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TABLE I: Computational complexity per pixel and iteration 


Algorithm 

IPSpLRU 

ADSpLRU 

CSUnSAl js 


MMV-ADMM |10| 

BilCE 


Computational complexity 

0{KN'^ + K^) 

0{KN‘^ + KLN) 

0{N^) 

0{KN^ + KLN) 

0{N^) 


IV. Experimental Results 

This section unravels the performanee charaeteristies of the proposed IPSpERU and AD- 
SpERU algorithms via experiments eondueted both on simulated and real data. We eompare our 
teehniques with three well-known state-of-the-art unmixing algorithms, namely, the nonnegative 
eonstraint sparse unmixing by variable splitting and augmented Eagrangian algorithm (CSUn- 
SAE), Q, the reeently reported nonnegative eonstraint joint-sparse method (MMV-ADMM), 
0 . and finally the fast Bayesian inferenee iterative eonditional expectations (BilCE) unmixing 
algorithm, [j^. The computational complexity (in terms of the number of multiplieations) of all 
the tested algorithms is given in Table |I| As shown in the Table, the spatial eorrelation-aware 
algorithms namely, IPSpERU, ADSpERU and MMV-ADMM, present higher complexity sinee 
the information from K pixels is used for the unmixng of a single pixel. Moreover, it is noticed 
that among the two proposed algorithms, IPSpERU has lower eomputational eomplexity than 
ADSpERU per iteration, resulting from its more simplistic incremental approach. 

In what follows, we first refer to the parameters’ setting established for all the involved algo¬ 
rithms, and the performanee evaluation metries that are utilized in the experimental proeedure. 
To eorroborate the effectiveness and robustness of the proposed algorithms we execute five 
different types of synthetic data experiments whose detailed description is given below. Einally, 
we empirieally eompare the abundanee maps as revealed by all examined algorithms, when 
applied on a real hyperspectral image. 

A. Setting of Parameters and Performance Evaluation Criteria 

Eor simplicity reasons, we use 7 for the sparsity imposing parameter in all tested algorithms 
(exeept BilCE whieh has no regularization parameters, 0), /i for the Eagrange multiplier 
regularization parameter of the ADMM-type teehniques and A for the relevant to /r regular¬ 
ization parameter of IPSpERU. Additionally, the low-rank promoting parameter of the proposed 
algorithms, is denoted by r. Parameters r and 7 are fine tuned with 10 different values, as 
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TABLE II: Parameters Setting 


Algorithm 

T (rank regularization parameter) 

7 (sparsity regularization parameter) 

M 

A 

IPSpLRU 

1 

o 

05 

1 

O 

o 

1 

O 

o' 

{0,10-1°, 10“®,. 

..,10-1} 

Not applicable 

0.5 

ADSpLRU 

{0, 

{0,10-1°, 10-°,. 

..,10-1} 

10-° 

Not applicable 

CSUnSAL 

Not applicable 

{0,10-1°, 10-°,. 

..,10-1} 

10-° 

Not applicable 

MMV-ADMM 

Not applicable 

{0,10-1°, 10-°,. 

..,10-1} 

10-° 

Not applicable 


shown in Table On the other hand, the Lagrange multiplier regularization parameter /i and 
the regularization parameter A of IPSpLRU, whieh influenee to a less extend the effieieney of 
the eorresponding algorithms, are set to a fixed value. In order to assess the performance of the 
proposed algorithms and the competing ones, for the experiments conducted on synthetic data 
we consider two metrics. First, the root mean square error (RMSE), 


RMSE = 


\ 


Nn 


E 

2=1 


|Wi - Wi||2, 


(30) 


where Wj and Wj represent the estimated and actual abundance vectors of the f-th pixel respec¬ 
tively, n is the total number of the pixels in the image under study, and N, as mentioned in 
previous sections, stands for the number of endmembers. The second metric, is the signal-to- 
reconstruction error (SRE), Q, which reflects the ratio between the power of the signal and the 
power of the estimation error, and is given by the following formula 

^Er=iiiw.iii 


SRE = lOlog 


10 


'.ill 


(31) 


iEr=illw*-Wi 

Of great importance is to notice that for the sliding window-based algorithms, the abundance vec¬ 
tors Wj’s and their estimates Wj’s coincide with the central column vectors of the corresponding 
abundance matrices Wj’s and their estimates W/s, as becomes clear from Fig. ^ 


B. Experiments on Simulated Datacubes 

In the sequel, N endmembers are randomly selected from the USGS library Q G 
p^ , so as to form our endmembers’ dictionary $. Their reflectance values correspond to L = 
224 spectral bands, uniformly distributed in the interval 0.4 — 2.5pm. The LMM of eq. (1) is 
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then utilized for generating speetral signatures subjeet to given, different in eaeh experiment, 
abundanee matriees W’s. 

1) Reweighting coefficients efficiency and convergence behavior of IPSpLRU and ADSpLRU: 
Herein, we aspire to demonstrate the merits emerging from the utilization of reweighting of A 


and b from (29), on the estimation performanee of the proposed algorithms. In light of this, we 
eonsider a rank 3 and sparsity level 10% (i.e., 10% of its entries are nonzero) abundanee matrix 
eorresponding to A = 50 endmembers, K = 9 pixels. Then, K = 9 speetral signatures are 
generated aeeording to the LMM and eontaminated by Gaussian noise sueh that SNR = 30dB. 

For p = 100 realizations. Fig. depiets the normalized mean squared estimation error 


-, where W* and Wj are the estimated at 


(NMSE) (defined as NMSE(f) = f Y7i=i 
the f-th iteration and the true matriees respeetively, of the z-th realization) as it evolves over 
2000 iterations. Three different eases are investigated, eorresponding to: a) updating weighting 
eoeffieients from ( [29] ) b) keeping fixed the weighting eoeffieients based on (28) and e) no 
weighting eoeffieients i.e. the weighted norms degenerate to their non-weighted versions by 
setting b = 1 and A = [1,1,...!]. As it is elearly evident in Fig. 1^ both IPSpLRU and 
ADSpLRU aehieve remarkably higher estimation aeeuraey in terms of NMSE, when using 
reweighting as eompared to the ease that fixed or no weights are employed. It is thus empirieally 
verified that the enhaneed effieieney of the reweighted ii and nuelear norms, emphatieally 


advoeated in [25|-[27|, is retained when using the sum of these two norms. The priee to be 
paid is that sueh an option might inerease numerieal risks, sinee the problem is rendered non- 
eonvex and (yet) no theoretieal eonvergenee analysis has been established. Nevertheless, it is 
worthy to mention that, despite the faet that eonvergenee is not theoretieally guaranteed, in all our 
experiments both IPSpLRU and ADSpLRU exhibited a very robust behavior in their eonvergenee 
proeess. 

It is also notieed that ADSpLRU needs less iterations to eonverge as eompared to IPSpLRU 
and it eonverges to a slightly lower NMSE. This results from the inherent nature of the two 
proposed algorithms, as explained before. Interestingly, the faster eonvergenee rate of ADSpERU 
with reweighting eomes at the priee of its higher per iteration eomputational eomplexity as 
eompared to that of IPSpERU. 

2) A toy example: In this experiment our goal is to highlight the signifieanee of the approaeh 
followed in this work, i.e., the simultaneous ineorporation of sparsity and low-rankness on the 
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Fig. 2: Convergence curves of IPSpLRU and ADSpLRU for a) updating weighting coefficients 
b) fixed weighting coefficients and c) no weighting coefficients. 


abundance estimation problem. To this end, we initially derive the single prior counterparts of 
our algorithms. We first focus on the low-rankness assumption, thus the sparsity imposing norm 
is ignored (7 = 0 ). IPSpLRU and ADSpLRU are then reduced to their modified versions, namely, 
IPLRU and ADLRU respectively. As implied by their names, the aforementioned methods allow 
exclusively for the low-rank assumption. Similarly, IPSpU and ADSpU are formed by accounting 
solely for sparsity. That said, IPSpU and ADSpU emerge after dropping the low-rank prior 
constraint (r = 0). Next, we generate a N x K (where = 50 and K = 9) simultaneously 
sparse and low-rank abundance matrix W of rank 2 with sparsity level 20%, which is graphically 


illustrated in Fig. 3a Using this W we generate the Lx K observations matrix Y via the LMM 
in Eq. ([T]), where the noise matrix E is Gaussian i.i.d. and SNR=35dB. 

Fig. shows the merits of the proposed IPSpLRU and ADSpLRU algorithms. Specifically, 
it appears that the concurrent exploitation of sparsity and low-rankness leads to significantly 
more accurate abundance matrix estimates, as compared to their single constraint counterparts. 
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2 4 6 8 
PIXELS 


(a) W, Ground truth 


Algorithm 

RMSE 

IPSpLRU 

0.0056 

IPSpU (sparse only) 

0.0121 

IPLRU (low-rank only) 

0.0098 

ADSpLRU 

0.0061 

ADSpU (sparse only) 

0.0130 

ADLRU (low-rank only) 

0.0102 



2 - 4-6 8 
PIXELS 



0.4 

0.3 

0.2 




2 4 6 8 
PIXELS 


0.12 


0.1 


0.08 


0.06 


0.04 


0.02 


O 
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PIXELS 
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0.3 


0.2 




2 4 6 8 
PIXELS 


I 

I 


0.1 2 


0.1 


0.08 


0.06 


0.04 


0.02 


O 


(b) W, IPSpLRU (c) residual, IPSpLRU 


(d) W, ADSpLRU (e) residual, ADSpLRU 
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I 


0.12 
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I 


0.08 


0.06 
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0.02 


O 


(f) W, IPLRU (g) residual, IPLRU 


(h) W, ADLRU (i) residual, ADLRU 
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PIXELS 




2 4 6 8 
PIXELS 


I 


0.1 2 


0.1 


I 


0.08 


0.06 


0.04 


0.02 


O 


(j) W, IPSpU 


(k) residual, IPSpU 


(1) W, ADSpU (m) residual, ADSpU 


Fig. 3: Sparse and low-rank algorithms versus their sparse only and low-rank only eounterparts. 


namely, IPLRU, IPSpU and ADLRU, ADSpU respectively. This is clearly seen in terms of the 
RMSE, as well as from a careful visual inspection of both the recovered abundance matrices 
and their residuals with the true abundace matrix (i.e. |W — W|), depicted in pair from Fig 


3b 


- Fig. 3m 


3) The key role of the parameters 7, r.- As explained earlier, parameters 7 and r control 
the imposition of sparsity and low-rankness, respectively, on the abundance matrix W. Herein, 


October 15, 2015 


DRAFT 



















































21 



(a) sparsity level=100%, (b) sparsity level=10%, (c) sparsity level=20%, (d) sparsity level=10%, 


rank=l, IPSpLRU 


rank=4, IPSpLRU 


rank=5, IPSpLRU 


rank=9, IPSpLRU 



(e) sparsity level=100%, (f) sparsity level=10%, (g) sparsity level=20%, (h) sparsity Ievel=10%, 


rank=l, ADSpLRU 


rank=4, ADSpLRU 


rank=5, ADSpLRU 


rank=9, ADSpLRU 


Fig. 4: RMSE as a function of the low-rankness and the sparsity regularization parameters r 
and 7 , respeetively. 


we unveil the dependeney of the optimal (with respeet to RMSE minimization) set of these 
parameters on the inherent strueture of the sought abundanee matrix. In this vein, five different 
types of abundanee matriees are generated, eaeh refleeting a speeifie eombination of rank and 
sparsity level. Next, K = 9 linearly mixed pixels are produeed, eorrupted with Gaussian i.i.d. 
noise and SNR=35dB. A number of 100 independent realizations is run for eaeh of the five 
experiments, and the average RMSE is demonstrated as a funetion of r and 7 . As shown in 
Eig. Q in the first ease (Eigs. 4a and 4e), whieh eorresponds to solely low-rank abundanee 
matriees (without any presenee of sparsity), the sparsity promoting parameter 7 does not affeet 


the estimation aeeuraey. In a similar manner, in the fourth experiment (Eigs. 4d and 4h), where 
the abundanee matrix is eonsidered full-rank and sparse, the low-rank promoting parameter has 
no impaet on the estimation performanee. Notably, in the other two oases (oolumns 2 and 3) 
where both sparse and low-rank abundanee matriees are eonsidered, RMSE is minimized for 
non-zero values of both r and 7 . Sueh a result is eonsistent with the fundamental premise of 
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SRE vs SNR yg gMR 




Fig. 5: Robustness to white noise (SRE & RMSE). 


our algorithms, whieh is the improvement in the abundance matrix estimation by simultaneously 
exploiting sparsity and low-rankness. 

Moreover, the above results indicate that the optimal choice of r, 7 depends on the particular 
structure (sparse and/or low-rank) of the abundance matrix. Thus, a proper selection of these pa¬ 
rameters shall involve fine-tuning schemes, which are commonplace when it comes to algorithms 
dealing with regularized inverse problems. 

4) Robustness to noise: In this experiment we aim at exhibiting the robustness of the proposed 
algorithms to white and correlated noise corruption. To this end, we stick with a specific 
simultaneously sparse and low-rank abundance matrix W of sparsity level 20% and rank 3. 
Based on this W, K = 9 linearly mixed pixels are generated, in the same way as described 
above. Then, depending on the case, white or colored Gaussian noise contaminates the data. 16 
SNR values are considered ranging from 10 to 40 dB, while 100 realizations are run for each 
SNR value, and the mean of the RMSE and SRE metrics is calculated. 

• White Gaussian Noise: Eig. shows the RMSE and SRE curves obtained for the proposed 
IPSpERU, ADSpERU and the three competing algorithms, namely, CSUnSAE, MMV- 
ADMM and BilCE. It is easily seen that both IPSpERU and ADSpERU attain remarkably 
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SRE vs SNR RMSE vs SNR 




Fig. 6: Robustness to colored noise (SRE & RMSE). 


better results comparing to CSUnSAE, MMV-ADMM and BilCE in all the examined 
SNR values. Additionally, we note that ADSpERU performs slightly better as compared 
to IPSpERU, especially for SNR values greater than 32dB. The price to be paid is that the 
computational complexity per iteration of ADSpERU is higher than that of IPSpERU. It is 
hence shown that sparse and low-rank methods are robust to different levels of white noise. 
At the same time, IPSpERU and ADSpERU outperform the sparse only CSUnSAE and 
BilCE algorithms and the joint-sparse MMV-ADMM algorithm, provided that both sparsity 
and low-rankness characterizes the abundance matrix. 

• Colored Gaussian Noise: Actually, in real hyperspectral images the noise that corrupts the 
data is rather structured than white. Thus, to assess the behavior of the proposed methods 
in such realistic conditions, we simulate correlated Gaussian noise that adds up to the 
linearly mixed pixels. Eig. [^illustrates the effectiveness of the tested algorithms in terms of 
RMSE and SRE, for different SNR values. Therein as well, we can see that IPSpERU and 
ADSpERU achieve better results than their competing algorithms in the whole range of the 
examined SNRs. Eurthermore, ADSpERU performs better for high SNR values (> 32dB), 
as compared to IPSpERU. As a result, the robustness of our proposed methods is also 
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row 

column 


2nd 

3rd 

4th 

joint sparse - 

(4.1) 

(8.2) 

(12,3) 

(16.4) 

low-rank - 2"'^ 

(100,1) 

(100.2) 

(100,3) 

(100,4) 

sparse & low-rank - 3'''* 

(4.2) 

(8.2) 

(12,2) 

(16,2) 

sparse & low-rank - 4*^ 

(4.3) 

(8.3) 

(12.3) 

(16.3) 


(b) Structure of W in each block of the synthetic image, each cell 
contains the pair : (sparsity-level%, rank(W).) 



Algorithm 

1®* row 

2"'' row 

3’'"' row 

4*^ row 

RMSE 

SRE 

RMSE 

SRE 

RMSE 

SRE 

RMSE 

SRE 

ADSpLRU 

0.009 

28.96 

0.078 

16.62 

0.032 

18.71 

0.029 

19.62 

IPSpLRU 

0.008 

28.39 

0.081 

16.41 

0.026 

21.01 

0.030 

19.81 

CSunSAL 

0.026 

19.81 

0.117 

12.39 

0.052 

13.88 

0.047 

14.99 

MMV-ADMM 

0.030 

18.00 

0.105 

12.99 

0.061 

12.32 

0.056 

13.16 

BilCE 

0.028 

21.71 

0.263 

6.72 

0.043 

17.83 

0.060 

15.81 


(c) RMSE (10 and SRE (dB) results on synthetic image for 
(a) Synthetic Image, 16 blocks of size 10 x 10 pixels each. s^ch row. 


Fig. 7: Structure of the synthetic image and results. 


corroborated in the presence of correlated noise with different magnitude. 

5) Synthetic Image: This experiment highlights the effectiveness of the proposed methods in 
estimating sparse, low-rank or both sparse and low-rank abundance matrices. Focused on this 
purpose we form a simulated hyperspectral image using the linear mixing model (1) and the same 
above-mentioned endmembers’ dictionary $. As shown in Fig. the simulated hyperspectral 
image consists of 4 rows each consisting of 4 10 x 10 blocks of pixels. Each of the “block rows” 
is generated by abundance matrices of a distinct structure. To be more specific, the first row 
is generated by joint-sparse W’s, the second by solely low-rank W’s, while rows 3 and 4 are 
produced by simultaneously sparse and low-rank abundance matrices. The pixels in each block 
correspond to abundance matrices of a particular combination of sparsity level and rank. The 


detailed description of these structures is depicted in the table of Fig. 7b The linearly mixed 


pixels are corrupted by white Gaussian i.i.d. noise such that SNR = 30dB. The table in Fig. 7c 


contains the obtained RMSE and SRE for all algorithms tested. It is worth pointing out that our 
introduced IPSpERU and ADSpERU algorithms outperform their rivals, not only in the “both 
sparse and low-rank” rows 3 and 4, but also in rows 1 and 2 that correspond to either sparse 
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(a) 5th PC of the Salinas Valley scene. 


1. Grapes 

2. Broccoli A 

3. Broccoli B 

4. Lettuce A 

5. Lettuce B 

6. Lettuce C 

7. Lettuce D 

8. Corn 


(h) Rough ground truth information for a part of the Salinas 
valley scene under study. 




(c) Spectral signatures of the 17 endmembers manually selected 
as pure pixels. 

Fig. 8: Salinas valley image and endmembers’ dietionary 


only or low-rank only W’s. 


C. Experiment on Real Data 

This section illustrates the performance of the proposed algorithms when applied on a real 
hyperspectral image. The hyperspectral scene under examination is a portion of the widely 
used Salinas vegetation scene acquired by AVIRIS sensor over Salinas Valley in California. 
This scene contains eight different vegetation species, namely grapes, brocolli_A, brocolli_B, 


lettuce_a, lettuce_b, lettuce_c, lettuce_d, com, as shown in Fig. 8b Salinas hyperspectral image 
consists of L = 204 spectral bands and its spatial resolution is 3.7 meters. Taking the principal 
components (PCs) of the image, it can be seen that only the first 6 of them contain meaningful 
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grapes 


brocolli_a 


brocolli_b 


corn 



(a) IPSpLRU 



(b) ADSpLRU 



(c) CSUnSAL 



(d) MMV-ADMM 



(e) BilCE 


Fig. 9: Abundance maps of Salinas hyperspectral image. 
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information. Focusing on them, we can see that the first PCs give more rough information about 
the formation of the vegetation, while less signifieant PCs give more refined information about 
the vegetation formation, |371. Fig. [8a| shows the 5th principal component (PC) of the scene under 
study, where most of the vegetation is depieted. The endmembers dictionary $ is composed of 
17 pure pixels’ spectral signatures, whieh have been seleeted manually, as in and depieted 


in Fig. 8e 


Fig. 1^ shows the abundanee maps corresponding to the region of interest, as obtained by 
the proposed IPSpLRU, ADSpLRU and the three state-of-the-art eompeting algorithms namely 
CSUnSAL and MMV-ADMM and BilCE for 7 = IQ-^, r = 10-^ A = 0.5 and /i = IQ-^. Four 
different maps are depieted, corresponding to four vegetation species, namely: grapes, broeolli a, 
broeolli b and eom. It is worth pointing out that sinee detailed ground truth information is 
not available, the evaluation is done in qualitative terms. From a careful visual inspeetion of 
the generated maps, we ean see that the abundanees obtained by IPSpLRU and ADSpLRU 
present patterns whieh are eloser to those revealed by the first five principal components of the 
hyperspeetral image provided in p8| . This is particularly clear for the maps corresponding to 
broeolli a and broeolli b. More specifically, it is shown that the presenee of these two species, 
whieh is mainly loeated in two distinet regions, is better emphasized by the proposed algorithms. 
In addition, the erroneous detection of these vegetation types is eliminated more effectively by 
IPSpLRU and ADSpLRU, as also verified by Figs 8a and Hence, it is eorroborated that the 
exploitation of the inherent spatial eorrelation existing in hyperspeetral images, ean lead us to 
qualitatively better results, thus verifying the signifieanee of our approaeh. 


V. Conclusions and Future Directions 

In this paper we presented a novel approaeh for performing hyperspeetral image unmixing 
exploiting simultaneously sparsity and spatial eorrelation. A novel eost function was first in¬ 
troduced eomprising a least squares proximity eomponent regularized by a linear eombination 
of the weighted £1 norm and the weighted nuelear norm of the latent abundance matrix. The 
unmixing problem was thus treated as a sparse redueed-rank regression problem. Two different 
algorithms were then developed for solving it, namely an incremental proximal type algorithm 
called IPSpLRU, and an ADMM based strategy ealled ADSpLRU. Extensive simulations on both 
synthetic and real data eorroborate the effeetiveness of the proposed approaeh and algorithms. 
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compared to other related state-of-the-art unmixing schemes. The derivation of more computa¬ 
tionally efficient schemes alleviating the need for SVD is under current investigation. Another 
relevant future research direction is the exploitation of the specific structure or pattern of sparsity 
in the abundance matrices imposed implicitly by the low-rankness property, which could further 
improve estimation performance. This is also a topic of interest in the framework of a future 
work. 
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